# Simulation of Romer-Rosenthal Game vs. Nash Bargaining Solution (NBS)
# JBS Feb 22 2013

rm(list=ls())
library(sfsmisc)
library(grid)

# Number of issues
nissues<-150

# Establish true locations of two actors, sq is assumed to be 0
x1 <- round(runif(nissues, min=0, max=60))
x2 <- round(runif(nissues, min=40, max=100))
x2 <- ifelse(x1==x2,x2+1,x2) # Need to make sure that the two positions do not agree

# Create True Romer-Rosenthal Outcome, X2 is agenda setter
outtrue <- ifelse((2*x1)<x2,2*x1,x2)

# Create NBS Function
# s is sq, a is position of x1 and b is position of x2. We seek to maximize c*d (combined utility) by varying position x
nb <- function(s,a,b,x){
			c <- ((-1*(a-x))^2)-((-1*(a-s))^2)
			d <- ((-1*(b-x))^2)-((-1*(b-s))^2)
			c*d		
		}

## True NBS outcome
outnb <- rep(NA,nissues)

for (i in 1:nissues){
	res<-optimize(nb,lower=min(x1[i],x2[i]),upper=max(x1[i],x2[i]),a=x1[i],b=x2[i],s=0,maximum=TRUE)
	outnb[i] <- res$maximum	
}

#########################
#########################
# Monte Carlo Simulations
# Measure measurement error same for both actors
# Run 1000 sims for each level of measurement error

# set levels of sd to vary
sd <- c(2,4,8,12)
#sd<-8

est <- vector("list",length(sd))
prop <- rep(NA,length(sd))
outobsRR<-vector("list",1000)
outobsNBS<-vector("list",1000)

for (j in 1:length(sd)){

esta <- matrix(NA,1000,2)

for (k in 1:1000){
	# Create observed positions, outcomes, and sq (our measurements)
	x1obs <- round(x1+rnorm(nissues,0,sd[j]),0)
	x2obs <- round(x2+rnorm(nissues,0,sd[j]),0)
	outobs <- round(outtrue+rnorm(nissues,0,sd[j]),0)
	sqobs  <- round(rnorm(nissues,0,sd[j]),0)
	
	#rescale to place observed obs on 0-100 scale
	absmax <- max(apply(cbind(x1obs,x2obs,outobs,sqobs),1,max))
	absmin<- min(apply(cbind(x1obs,x2obs,outobs,sqobs),1,min))
	range<- absmax-absmin
	x1obs <- round(((x1obs+ abs(absmin))/range)*100,0)
	x2obs <- round(((x2obs+ abs(absmin))/range)*100,0)
	outobs<- round(((outobs+ abs(absmin))/range)*100,0)
	sqobs <- round(((sqobs+ abs(absmin))/range)*100,0)
	x2obs <- ifelse(x1obs==x2obs,x2obs+1,x2obs) # Need to make sure that the two positions do not agree
	
	# For simplicity, ensure that sq is never to left of both observed positions.
	sqobs[sqobs>x1obs & sqobs>x2obs] <- x2obs[sqobs>x1obs & sqobs>x2obs] - 1

	# Create theoretical outcomes based on observed positions
	# Create R-R prediction ---- SQ observed with error
			outrrobs <- rep(NA,nissues) # Holding bin for R-R observations
			# create scenarios --- SQ inbetween observed positions
			outrrobs[x1obs<=sqobs & sqobs<=x2obs | x2obs<=sqobs & sqobs<=x1obs]<-sqobs[x1obs<=sqobs & sqobs<=x2obs | x2obs<=sqobs & sqobs<=x1obs]
			# create scenarios --- SQ to the right of both positions
			outrrobs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)>=(x2obs-sqobs)]<-x2obs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)>=(x2obs-sqobs)]
			outrrobs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)<(x2obs-sqobs)]<-2*x1obs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)<(x2obs-sqobs)]
			# create scenarios --- SQ to the left of both positions
			
			outobsRR[[k]]<-outrrobs
			
	outnbobs <- rep(NA,nissues)
			
	for (i in 1:nissues){
		res<-optimize(nb,lower=min(x1obs[i],x2obs[i]),upper=max(x1obs[i],x2obs[i]),a=x1obs[i],b=x2obs[i],s=sqobs[i],maximum=TRUE)
		outnbobs[i] <- res$maximum	
	}
				outobsNBS[[k]]<-outnbobs
	# Run model with theoretical outcomes predicting actual outcomes.
	mod1 <- lm(outobs ~ outrrobs + outnbobs - 1)
	esta[k,]<-mod1$coefficients
}
	est[[j]]<-esta
	prop[j]<-mean(ifelse(esta[,1]>esta[,2],1,0))
}

allest<-as.matrix(cbind(est[[1]],est[[2]],est[[3]],est[[4]]))

colnames(allest)<-c("RR","NB","RR","NB","RR","NB","RR","NB")
boxplot.matrix(allest,ylab="Coefficients",ylim=c(-0.3,1.50))
text(x = 1.4, y = 1.3, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=1.4, y=1.4, "Model 1 SD2")
text(x = 3.4, y = 1.3, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=3.4, y=1.4, "Model 2 SD4")
text(x = 5.4, y = 1.3, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=5.4, y=1.4, "Model 3 SD8")
text(x = 7.4, y = 1.3, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=7.4, y=1.4, "Model 4 SD12")


m1percor<-mean(ifelse(allest[,1]>allest[,2],1,0))
m1percor
m2percor<-mean(ifelse(allest[,3]>allest[,4],1,0))
m2percor
m3percor<-mean(ifelse(allest[,5]>allest[,6],1,0))
m3percor
m4percor<-mean(ifelse(allest[,7]>allest[,8],1,0))
m4percor


merrorRR<-rep(NA,1000)
merrorNB<-rep(NA,1000)

for (i in 1:1000){
	merrorRR[i]<-mean(abs(outtrue-outobsRR[[1]]))
	merrorNB[i]<-mean(abs(outnb-outobsNBS[[1]]))
}

mean(merrorRR)
mean(merrorNB)
########################################
########################################
# MC Simulation RR vs Average Position


sd <- c(2,4,8,12)

est <- vector("list",length(sd))
prop <- rep(NA,length(sd))

for (j in 1:length(sd)){

esta <- matrix(NA,1000,2)

for (k in 1:1000){
	# Create observed positions, outcomes, and sq (our measurements)
	x1obs <- x1+rnorm(nissues,0,sd[j])
	x2obs <- x2+rnorm(nissues,0,sd[j])
	outobs <- outtrue+rnorm(nissues,0,sd[j])
	sqobs  <- rnorm(nissues,0,sd[j])
	
	#rescale to place observed obs on 0-100 scale
	absmax <- max(apply(cbind(x1obs,x2obs,outobs,sqobs),1,max))
	absmin<- min(apply(cbind(x1obs,x2obs,outobs,sqobs),1,min))
	range<- absmax-absmin
	x1obs <- round(((x1obs+ abs(absmin))/range)*100,0)
	x2obs <- round(((x2obs+ abs(absmin))/range)*100,0)
	outobs<- round(((outobs+ abs(absmin))/range)*100,0)
	sqobs <- round(((sqobs+ abs(absmin))/range)*100,0)
	x2obs <- ifelse(x1obs==x2obs,x2obs+1,x2obs) # Need to make sure that the two positions do not agree
	 
	# Ensure that sq is never to left of both observed positions.
	sqobs[sqobs>x1obs & sqobs>x2obs] <- x2obs[sqobs>x1obs & sqobs>x2obs] - 1

	# Create theoretical outcomes based on observed positions
	# Create R-R prediction ---- SQ observed with error
			outrrobs <- rep(NA,nissues) # Holding bin for R-R observations
			# create scenarios --- SQ inbetween observed positions
			outrrobs[x1obs<sqobs & sqobs<x2obs | x2obs<sqobs & sqobs<x1obs]<-sqobs[x1obs<sqobs & sqobs<x2obs | x2obs<sqobs & sqobs<x1obs]
			# create scenarios --- SQ to the right of both positions
			outrrobs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)>(x2obs-sqobs)]<-x2obs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)>(x2obs-sqobs)]
			outrrobs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)<(x2obs-sqobs)]<-2*x1obs[x2obs>sqobs & x1obs>sqobs & 2*(x1obs-sqobs)<(x2obs-sqobs)]
			# create scenarios --- SQ to the left of both positions

	# This version just uses the average position rather than NBS
	outnbobs <- (x1obs+x2obs)/2

	# Run model with theoretical outcomes predicting actual outcomes.
	mod1 <- lm(outobs ~ outrrobs + outnbobs - 1)
	esta[k,]<-mod1$coefficients
}
	est[[j]]<-esta
	prop[j]<-mean(ifelse(esta[,1]>esta[,2],1,0))
}

allest<-as.matrix(cbind(est[[1]],est[[2]],est[[3]],est[[4]]))

colnames(allest)<-c("RR","AVE","RR","AVE","RR","AVE","RR","AVE")
boxplot.matrix(allest,ylab="Coefficients",ylim=c(-0.3,2))
text(x = 1.4, y = 1.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=1.4, y=1.9, "Model 1 SD2")
text(x = 3.4, y = 1.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=3.4, y=1.9, "Model 2 SD4")
text(x = 5.4, y = 1.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=5.4, y=1.9, "Model 3 SD8")
text(x = 7.4, y = 1.8, '{', srt = 270, cex = 8, family = 'Helvetica Neue UltraLight')
text(x=7.4, y=1.9, "Model 4 SD12")


m1percor<-mean(ifelse(allest[,1]>allest[,2],1,0))
m1percor
m2percor<-mean(ifelse(allest[,3]>allest[,4],1,0))
m2percor
m3percor<-mean(ifelse(allest[,5]>allest[,6],1,0))
m3percor
m4percor<-mean(ifelse(allest[,7]>allest[,8],1,0))
m4percor
############################################
